function wave_movie
%
%  plots sine series solution of the wave equation
%       diff(u,x,x) = diff(u,t,t)   for 0 < x < L, 0 < t < tmax
%  where
%      u = 0  at x=xL,xR  and  u = g(x)  and  u_t = h(x)  at t = 0

tmax=6;
L=1;

% the code begins
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

v = VideoWriter('/Users/mark/Desktop/wave','MPEG-4');
open(v)

% get(gcf)
set(gcf,'Position', [1 996 692 349])
subaxis(1,1,1,1,'MT',0.08,'MB',0.1,'MR',-0.01,'ML',0.1,'P',0.04)

N=40;
M=150;
frames=3;

% generate the points along the x-axis, x(1)=xL and x(N)=xR
x=linspace(0,L,N);

% generate the points along the t-axis, t(1)=0 and t(M)=tmax
t=linspace(0,tmax,M);


s=zeros(N,M);
for j=1:M
    for ix=1:N
        sum=0;
        for n=1:50
            sum=sum+(a(n)*cos(n*pi*t(j)/L)+b(n)*sin(n*pi*t(j)/L))*sin(n*pi*x(ix)/L);
        end
        s(ix,j)=sum;
        
    end
end

yM=max(max(s))
ym=min(min(s))

for j=1:M
    clf
    subaxis(1,1,1,1,'MT',0.08,'MB',0.1,'MR',-0.01,'ML',0.06,'P',0.04)
    axis([0 L ym yM])
    hold on
    
    % say=['Solution of Wave Equation: \quad $u(x,0) = x^2(1-x)^2$ \quad $u_t(x,0) = 0$'];
    %say=['Solution of Wave Equation: \quad $u(x,0) = x^2(1-x)^2(x-0.25)$ \quad $u_t(x,0) = 0$'];
    say=['Solution of Wave Equation: \quad $u(x,0) = x^2(1-x)^2\sin(17x/3)$ \quad $u_t(x,0) = 0$'];
    title(say,'FontSize',14,'FontWeight','bold','Interpreter','Latex')
    
    xlabel('x-axis')
    ylabel('Solution')
    
    box on
    
    plot(x,s(:,j),'-b','LineWidth',1.6)
    if j>1
        plot(x,s(:,1),'--b','LineWidth',1)
    end
        
    say=['t = ',num2str(t(j),'%4.1f')];
    text(0.85*L,0.8*yM,say,'FontSize',20,'FontWeight','bold')
    
    set(gca,'FontSize',15,'FontWeight','bold')
    
    % make movie frame
    F = getframe(gcf);
    for i=1:frames
        writeVideo(v,F);
    end
    
    hold off
    
end

function y=a(n)
%y=4*(n^2*pi^2*(-1)^n - n^2*pi^2 - 12*(-1)^n + 12)/(n^5*pi^5);      % g=x^2(1-x)^2
%y=(3*n^2*pi^2*(-1)^n + n^2*pi^2 - 132*(-1)^n - 108)/(n^5*pi^5);      % g=x^2(1-x)^2(x-0.25)
y=972*n*pi*(243*pi^6*(-1)^n*sin(17/3)*n^6 - 33048*pi^4*(-1)^n*cos(17/3)*n^4 ...
+ 4887*pi^4*(-1)^n*sin(17/3)*n^4 - 33048*n^4*pi^4 - 2189175*pi^2*(-1)^n*sin(17/3)*n^2 ...
+ 34076568*(-1)^n*cos(17/3) + 9103789*(-1)^n*sin(17/3) + 34076568)/((3*pi*n + 17)^5*(3*pi*n - 17)^5); % g=x^2(1-x)^2(x-0.25)sin(17x/3)

function y=b(n)
y=0;


















